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Abstract 

Until  recently,  structure  and  properties  of  molecular 
crystals  could  not  be  predicted  computationally  since  the 
forces  acting  between  the  constituent  molecules  in 
crystals  were  not  known  sufficiently  accurately.  This 
situation  has  changed  with  the  development  of  an 
electronic  structure  method  called  symmetry-adapted 
perturbation  theory  based  on  the  density-functional 
description  of  monomers  [SAPT(DFT)] .  This  method  is 
sufficiently  efficient  to  be  applied  to  interactions  of 
energetic  molecules,  for  example  to  cyclotrimethylene 
trinitramine  (RDX).  Systems  even  larger  than  RDX  can 
be  treated,  for  example,  interaction  energies  for  the  dimer 
of  perylene,  containing  64  atoms,  have  been  computed 
using  SAPT(DFT).  The  SAPT(DFT)  potential  developed 
for  RDX  was  used  in  a  combined  molecular  packing, 
lattice  energy’  minimization,  and  molecular  dynamics 
approach  to  characterize  low-energy’  polymorphs  of  the 
RDX  crystal.  The  lowest-energy’  structure  corresponded 
to  the  observed  crystal  and  the  results  obtained  for  high- 
density  polymorphs  provide  new  information  on  the 
polymorphism  of  RDX.  The  SAPT(DFT)  method  should 
find  important  applications  in  development  of  new 
energetic  materials,  including  crystal  design,  screening 
molecules  for  co-crystallization,  and  identification  of  low- 
energy’  polymorphs. 

1.  Introduction 

Energetic  material  and  many  other  compounds  of 
interest  to  national  defense  are  crystals  of  organic 
molecules,  bound  by  intermolecular  (van  der  Waals) 
forces  and  dominated  by  the  dispersion  component  of 
these  forces.  The  knowledge  of  these  force  fields  is 


sufficient  to  predict  properties  of  such  crystals.  Although 
intermolecular  interaction  energies  can  be  computed  using 
electronic  structure  methods,  such  methods  are  currently 
too  costly  to  be  applied  to  compounds  of  interest  in  the 
field  of  energetic  materials.  Therefore,  the  predictions  for 
such  systems  had  to  be  based  on  empirical  potentials 
fitted  to  available  experimental  data  on  molecular 
crystals.  This  limits  such  predictions  to  systems  used  for 
the  parametrization  (in  fact,  also  to  near-experimental 
thermodynamic  conditions).  In  recent  years,  the 
predictive  power  of  theory  has  been  examined  in  blind 
tests  conducted  by  the  Cambridge  Crystallographic  Data 
Center1 1  and  found  to  be  very  low.  The  situation  had  not 
improved  in  subsequent  tests,  in  fact,  the  success  rate  of 
the  third  test121  was  lower  than  that  of  the  previous  ones. 
Recently,  a  fourth,  yet  unpublished  test  was  performed, 
and  the  success  rate  has  apparently  become  somewhat 
better. [3,4]  This  is  an  unsatisfactory  situation  in  view  of 
the  importance  of  such  crystals.  For  example,  energetic 
materials  are  crystals  of  large  organic  molecules  and 
polymorphism  of  drugs  is  a  major  problem  in 
pharmaceutical  industry.  This  inability  of  theory  to 
predict  properties  of  molecular  crystals  has  been 
considered  to  be  one  of  the  major  issues  in  modern 
molecular  science.  [5~8] 

An  important  question  for  the  predictions  of  crystal 
structure  is  how  accurate  the  force  fields  should  be.  The 
failures  of  predictions  discussed  above  indicate  that  the 
accuracy  of  the  force  fields  used  was  not  sufficient. 
Highly  accurate  force  fields  can  be  computed  ab  initio 
using  wave-function  (WF)  based  methods,  but  only  for 
molecules  significantly  smaller  than  those  forming 
molecular  crystals.  Density  functional  theory  (DFT)  can 
be  applied  to  systems  containing  hundreds  of  atoms,  but 
conventional  DFT  methods  cannot  describe 
intermolecular  interactions  dominated  by  the  dispersion 
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component.  Recently,  a  method  has  been  developed 
which  combines  WF -based  symmetry-adapted 
perturbation  theory  (SAPT)19'101  of  intermolecular 
interactions  with  the  Kohn-Sham  DFT  representation  of 
monomers111-171.  This  method  is  denoted  as  SAPT(DFT) 
or  SAPT-DFT.  The  efficiency  of  the  original 
formulations  of  SAPT(DFT)  has  been  increased113'18-201  by 
applying  the  density-fitting  method1211.  [For  the  current 
capabilities  of  SAPT(DFT),  and  references  to  earlier 
work,  see  References  18  and  19.]  In  SAPT(DFT),  the 
electron  correlation  effects  inside  monomers  are 
accounted  for  using  the  DFT  approach  which  results  in  a 
very  low-cost  method:  SAPT(DFT)  calculations  take  less 
time  than  supermolecular  DFT  calculations  for 
interactions  of  monomers  containing  about  20  atoms.  At 
the  same  time,  interaction  energies  computed  using 
SAPT(DFT)  are  as  accurate  as  those  computed  using 
high-level  WF-based  methods. 118,19,221  As  an  example, 
SAPT(DFT)  has  been  applied  to  compute  the  complete 
potential  surface  of  the  water  dimer1231  and  of  the  benzene 
dimer1221,  in  each  case  giving  results  in  excellent 
agreement  with  experimental  data.  These  results  indicate 
that  the  force  fields  based  on  SAPT(DFT)  should  be 
accurate  enough  to  predict  crystal  structures.  SAPT(DFT) 
has  recently  been  applied  to  several  large  systems  such 
the  pyrene  dimer1241  and  the  perylene  dimmer1251. 

2.  RDX  CRYSTAL 

In  order  to  test  the  capabilities  of  SAPT(DFT),  we 
have  decided  to  investigate  a  crystal  built  of 
cyclotrimethylene  trinitramine  (RDX)  monomers  under 
essentially  the  same  conditions  as  in  the  blind  tests,  i.e., 
without  any  use  of  experimental  information  about  the 
crystal  except  for  the  monomer’s  geometry  within  the 
crystal  taken  from  measurements  of  Reference  26.  The 
ab  w/7/o-optimizcd  monomer  geometry  is  very  close  to 
the  geometry  observed  in  crystals1271,  so  that  the  use  of  the 
former  geometry  would  give  essentially  identical  results. 
One  should  point  out,  however,  that  there  exists  a 
monomer  structure  with  energy  very  close  to  that  of  the 
global  minimum  geometry.1271  This  structure  could  not  be 
excluded  a  priori  to  be  the  one  appearing  in  the  crystal 
and  in  a  completely  first-principles  approach  one  should 
develop  an  RDX  dimer  potential  also  for  this  monomer’s 
structure. 

A  simplified  version  of  the  SAPT(DFT)  RDX  dimer 
potential  of  Reference  28  has  been  used  to  generate  a 
large  number  of  RDX  crystal  polymorphs  using  the 
molecular  packing  and  lattice  minimization  methods 
available  in  the  MOLPAK/WMIN  computer  code 
packages.1291  The  development  of  a  simplified  version 
was  necessary  since  WMIN  codes  cannot  use  the  original 
form  of  the  potential.1281  The  simplified  version  was 


obtained  by  refitting  the  atom-atom  potentials  uab  with  the 
standard  Buckingham-type  formula: 
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where  rab  is  the  atom-atom  distance  and  qa,  qb  are  the 
atomic  charges  of  Reference  28.  The  sites  were  located  at 
the  atomic  positions  of  the  monomers.  The  parameters 
C£b  were  fitted  to  the  long-range  only  SAPT(DFT) 
interaction  energies  using  combining  rules: 
£ob  _  ^caagbb  parameters  ^  and  fiab  were  fitted 

to  the  complete  set  of  SAPT(DFT)  interaction  energies 
using  the  numerical  procedure  described  in  Reference  28. 
The  simplification  of  the  potential  had  only  a  small 
impact  on  crystal  predictions.  We  can  judge  this  based  on 
the  results  of  isothermal-isostress  (NsT)  molecular 
dynamics  (MD)  simulations  (see  below)  which  could  be 
performed  with  both  potentials:  the  cell  vectors  of  the 
time-averaged  geometry  differed  by  0.2-0. 7%  in  the  two 
cases. 

The  molecular  packing  part  of  the  MOLPAK/WMIN 
codes  sampled  51  most  common  space-symmetry  groups 
to  produce  50,653  candidate  structures  for  each  group. 
For  the  500  most  dense  structures,  the  lattice  parameters 
and  monomer  orientational  parameters  were  optimized 
(within  restrictions  of  a  given  space  group)  using  WMIN. 
From  these  500  polymorphic  structures,  the  13  lowest- 
energy  structures  have  been  selected  for  the  NsT  MD 
simulations  at  ambient  conditions.  The  NsT  MD 
simulations  were  performed  using  the  original 
SAPT(DFT)  RDX  dimer  potential.1281  No  symmetry 
constraints  were  used  in  this  step.  Note  that  this  potential 
has  already  been  used  in  Reference  28  in  NsT  MD 
simulations  for  the  crystal  of  RDX,  but  only  for  the 
experimentally  observed  polymorph,  i.e.,  no  crystal 
structure  predictions  were  attempted.  After  the 
equilibration  phase,  the  simulations  were  used  to  generate 
time-averaged  information  about  the  RDX  crystal 
structure.  This  was  accomplished  by  averaging  the 
molecular  parameters  for  each  of  the  symmetry- 
equivalent  molecules  in  all  unit  cells  within  the  simulation 
supercell.  For  comparison,  similar  calculations  were 
performed  with  the  empirical  RDX  dimer  potential  of 
Sorescu,  Rice,  and  Thompson  (SRT)[30]  (fitted  on  the 
properties  of  the  RDX  crystal)  and  the  11  lowest-lying 
configurations  were  selected  for  the  molecular  dynamics 
simulations. 

In  Figure  1,  a  plot  of  lattice  potential  energies  and 
densities  of  the  resulting  crystal  structures  is  presented. 
Our  lowest  energy  structure  is  essentially  identical  to  the 
experimental  structure.  It  is  separated  from  other 
structures  by  a  fairly  substantial  gap,  so  that  our 
prediction  is  unambiguous.  The  agreement  achieved  by 
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our  first-principle  approach  is  even  better  than  that 
achieved  by  the  empirical  SRT  potential  despite  the  fact 
that  the  experimental  information  we  compare  to  was 
used  in  fitting  of  this  potential.  Although  the  SRT 
potential  does  predict  the  correct  lowest-energy  structure, 
its  predictions  for  other  polymorphs  are  very  different 
from  the  SAPT(DFT)  predictions:  the  pattern  of  points  in 
Figure  1  is  quite  different  in  the  two  cases  and  the 
distributions  of  points  cannot  be  made  to  overlap  by 
translations.  Since  SAPT(DFT)  description  of  any  low- 
energy  polymorph  should  be  equally  accurate  as  the 
description  of  the  lowest-energy  polymorph,  we  conclude 
that  the  empirical  potentials  are  not  well  suited  to 
investigate  polymorphism  of  crystals.  The  main  reason  is 
that  the  empirical  potentials  may  describe  well  only  the 
regions  of  the  configuration  space  which  are  probed  by 
the  lowest  energy  crystal  for  which  experimental 
information  is  available.  As  shown  in  Reference  28,  this 
is  the  case  for  the  RDX  dimer:  the  SRT  and  SAPT(DFT) 
potentials  are  very  close  to  each  other  for  the  crystal 
closest  neighbor  dimer  configuration,  but  not  for  other 
configurations  corresponding  to  various  local  minima. 
Since  the  latter  configurations  dominate  the  interactions 
in  some  polymorphs,  the  SRT  potential  cannot  describe 
these  polymorphs  well. 
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Figure  1.  Correlation  between  the  density  and  lattice  energy 
per  molecule  for  the  13  polymorphic  structures  of  the  RDX 
crystal  from  MD  simulations  using  the  SAPT(DFT)  potential 
and  comparison  to  simulations  using  the  empirical  SRT 
potential  of  Reference  30.  The  point  corresponding  to  the 
observed  crystal  is  denoted  by  a  red  up-triangle. 

The  excellent  agreement  of  SAPT(DFT)  predictions 
for  the  RDX  crystal  with  experiment  indicates  that  this 
method  can  be  reliably  used  for  other  molecular  crystals. 
Such  predictions  can  be  done  entirely  from  first  principles 
so  that  the  method  can  be  used  even  for  compounds  for 
which  no  experimental  data  are  available  and  methods 
based  on  empirical  force  fields  cannot  be  used. 
Therefore,  this  method  should  find  applications  in  crystal 
design,  in  particular  in  development  of  novel  energetic 


materials  and  (opto)electronic  devices,  screening 
molecules  for  co-crystallization,  and  identification  of  low- 
energy  polymorphs  of  pharmaceutical  compounds. 


Figure  2.  Dimer  of  TAP  UBA  (chain  A)  and  nucleoporin 
peptide  (chain  B) 

3.  Linearly-Scaling  SAPT(DFT) 

Whereas  many  of  the  applications  mentioned  above 
can  be  pursued  with  the  current  version  of  the 
SAPT(DFT)  codes,  the  present  scaling  of  SAPT(DFT),  as 
the  fifth  power  of  system  size,  will  not  allow  applications 
to  dimers  containing  much  more  than  about  50  atoms. 
Flowever,  for  large  systems  the  scaling  can  be  reduced  to 
linear  if  calculations  are  properly  localized  to  only  some 
fragments  of  such  systems.  A  general  strategy  for  a  linear 
SAPT(DFT)  approach  was  outlined  in  Reference  31.  This 
paper  also  presented  applications  of  this  approach  to  the 
electrostatic  component  of  interaction  energies  for  cases 
with  a  localized  intermolecular  contact  region.  Recently, 
the  method  has  been  extended  to  dimers  with  extended 
contact  areas  such  as  two  parallel  polymers  or  two 
parallel  planar  molecules.  The  new  algorithm  groups 
atoms  of  the  monomers  into  domains  and  calculates 
interdomain  contributions  either  with  the  electron  density 
overlap  effects  or  from  simple  asymptotic  formulas.  This 
approach  should  enable  us  to  calculate  electrostatic 
interactions  of  monomers  containing  several  hundred 
atoms.  Using  this  algorithm,  we  have  reproduced  the 
exact  electrostatic  interaction  energies  of  dimers 
containing  two  long  parallel  alkane  chains,  e.g., 
C36FI74  C36FI74,  to  within  about  one  percent.  The 
approximate  method  was  then  applied  to  systems  for 
which  exact  calculations  are  not  possible.  One  example 
of  such  application,  currently  under  development,  is 
shown  in  Figure  2.  No  hydrogens  are  displayed  in  either 
panel.  The  left  panel  shows  all  heavy  atoms,  whereas,  the 
right  panel  only  the  atoms  of  the  backbone.  The  dimer 
shown  in  Figure  2  is  a  biomolecule  with  known 
crystallographic  structure. |32J  It  is  a  complex  between  the 
transporter  associated  with  antigen  processing  (TAP),  also 
known  as  nuclear  export  factor  (NFX1),  ubiquitin- 
associated  (UBA)  monomer  (chain  A)  and  FVFG 
nucleoporin  peptide  (chain  B),  where  F  stands  of 
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phenylalanine  and  FG  for  phenylalanine-glycine  peptides. 
These  biomolecules  are  important  in  the  field  of  nuclear 
transport.1'2  331 
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